%% Historical Decomposition

%Which Variables should be cumulated for the Historical Decomp 
cum=n_price;


[nstr, TT, ndraws]=size(strerr);        
histirf(:,:,:,:)=zeros(n_var,TT,nstr,tel);

J = zeros(n_var,n_var*p); J(:,1:n_var) = eye(n_var);

for tel3=1:tel
   AA=eye(size(A,1));
   start=1;
   constant=0;
 for tt=1:TT-start+1
     
     constant=constant+   J* AA*J'*gamma_g(n_var*p+1,:,1)'; %geometric sum of intercepts,  converges to C*inv(I-A) 
        AA=A(:,:,1)*AA;

for strshoc=1:nstr;

    helpirf=squeeze(irfstr(1:n_var,strshoc,:,tel3)) ;
    helpirf2=[zeros(n_var,tt-1),   helpirf(:,1:TT-tt+1).*repmat(strerr(strshoc,tt+start-1,tel3),n_var,TT-tt+1) ];
    histirf(:,:,strshoc,tel3)= histirf(:,:,strshoc,tel3)+ helpirf2;
    histempt(:,tt,strshoc,tel3)=zeros(n_var,1)+constant;  %non shock componenent (initial conditions + constant)
     for ll=1:nlags;
     histempt(:,tt,strshoc,tel3)= histempt(:,tt,strshoc,tel3)+AA(1:n_var,(ll-1)*n_var+1:ll*n_var)*y(start+nlags-ll,:)' ;
     end
     
   end
   end
end


histirf_temp=histirf;
y_hist=y(start+nlags:end,:);
y_hist(:,cum )=cumsum(y(start+nlags:end,cum),1);
histirf(cum,:,:,:)=cumsum(histirf_temp(cum,:,:,:),2);
histempt(cum,:,:,:)=cumsum(histempt(cum,:,:,:),2);




%% Renormalize
%Set time axis on charts
startdate =size(varnames,1)-gapstart;                        
timestr = varnames(startdate+1:end,1); 
step=6*4;% Retrieve dates
ll=size(timestr,1);
begsp=(ll/step-floor(ll/step))*step;
spacing = begsp:step:size(timestr,1);                                    % Set end date on axes
spacing=floor(spacing);

IndexC = strfind(timestr,DateZero);
datelow = find(not(cellfun('isempty',IndexC)));
if size(histirf,2)~=size(timestr,1)
    err
end

IndexC = strfind(timestr,BoomEnd);
datehigh = find(not(cellfun('isempty',IndexC)));

IndexC = strfind(timestr,BustEnd);
datebust = find(not(cellfun('isempty',IndexC)));


for ihist=1:size(histirf,2)
histirf2(n_decomp,ihist,:,:)=histirf(n_decomp,ihist,:,:)-histirf(n_decomp,datelow,:,:);
histempt2(n_decomp,ihist,:,:)=histempt(n_decomp,ihist,:,:)-histempt(n_decomp,datelow,:,:);
end

y_hist(:,n_decomp)=y_hist(:,n_decomp)-ones(size(y_hist,1),1)*y_hist(datelow,n_decomp);
histirf(n_decomp,:,:,:)=histirf2(n_decomp,:,:,:);
histempt(n_decomp,:,:,:)=histempt2(n_decomp,:,:,:);
histirf_med=median(histirf(:,:,:,:),4);
histempt_med=median(histempt(:,:,:,:),4);
histirf_up=prctile(histirf(:,:,:,:),0.84*100,4);
histirf_down=prctile(histirf(:,:,:,:),0.16*100,4);
